Multi-scale Graph-based Guided Filter for De-noising Cryo-Electron Tomographic Data

Cryo-Electron Tomography is a leading imaging technique in structural biology, which is capable of acquiring two-dimensional projections of cellular structures at high resolution and close-to-native state. Due to the limited electron dose the resulting projections exhibit extremely low SNR and contrast. The 3D structure is then reconstructed and passed through a number of post-processing steps including de-noising and subtomogram averaging to provide a better understanding and interpretation. As CET is mainly used for imaging fine scale structures, any denoising method applied to CET images should be scale selective and in particular be able to preserve such fine scale structures. In this context, we propose a new denoising framework based on regularized graph spectral filtering with a full control of scale-space and global consistency. Using the gold-standard metrics, we show that our denoising algorithm significantly outperforms the state-of-the-art methods such as NAD, NLM and RGF in terms of noise removal and structure preservation. Cryo-electron tomography (CET) is a powerful imaging technique in biological sciences which bridges the gap between the molecular and the cellular structural biology [16], giving a better understanding of protein interactions and thus better drug delivery strategies. In principle, similar to Computed Tomography (CT) in Medical Imaging, CET acquires twodimensional projections at high resolution (around 20-50 Angstrom) of three dimensional (3D) cellular structures (called tomograms) at cryogenic (freezing) conditions under near-tonative state. Due to the low electron dose, necessary to avoid biological specimen damage, and limited tilt angle (typically ±60◦ to 70◦), a noisy (SNR typically 0.1 to 0.01) and extremely weak signal (low contrast) is formed in the resulting projections. These unfiltered projection images are then projected back to build the tomogram. In the reconstruction phase, the noise is propagated through the tomogram making the noise model more complicated. c © 2015. The copyright of this document resides with its authors. It may be distributed unchanged freely in print or electronic forms. 2 ALBARQOUNI et al.: MULTI-SCALE GRAPH-BASED GUIDED FILTER (MG2F) The reader is referred to [6] for more details on image formation and noise model in Electron Microscopy. Therefore, post-processing steps, such as noise reduction, after 3D reconstruction are necessary to provide a better visualisation and interpretation of the structure under scrutiny. However, this process is critical and could lead to wrong interpretation by erroneously removing fine structural information, that can not be discriminated from the noise. Conventional linear filters such as Gaussian kernels succeed in reducing noise, however, at the expense of blurring edges. Popular Non-linear anisotropic diffusion (NAD) [15] and its extended versions, which can be interpreted in terms of scale space theory, are extensively used in CET community due to their successful performance. However, NAD requires the diffusivity to be chosen carefully, which is sometimes quite challenging, and needs many iterations to converge. Non-local means (NLM) filter and its fast version[2] are investigated for denoising tomograms having redundant information, however, their performance is degraded when the window size is increased (spatial neighbourhood), in particular for high resolution CET data (> 20482). To date, both techniques are still used in the context of CET, however, advanced filtering algorithms which are able to smooth the noise while preserving the edges to increase the contrast as well are still in demand and highly desirable [11]. Bilateral filter (BF) [19] and its vectorized extension [14] have been successfully applied in computer vision community. One related technique is the rolling guidance filter (RGF) [5], which can be interpreted in terms of joint bilateral filter. However, it uses the filtered image as guidance rather than the original image which is commonly used in guided filters. This way, it succeeds in preserving the edges while smoothing the background. Another related paper to our work is [7], where bilateral filtering (BF) is interpreted in graph spectral domain addressing some open issues in [17] regarding emerging signal processing on graphs. As mentioned before, denoising CET images requires a proper scale selection as well as the preservation of fine scale structures. The proposed method is thus based on the following considerations: • By using a multi-scale pyramid for guidance we are able to detect meaningful scales and use them for guidance without oversmoothing fine scale structures. • Using a patch-based approach, we can take advantage of redundant structures in the whole image rather than using a pre-defined spatial window for averaging similar pixels or patches. This way, we can preserve the local and global consistencies. • By deriving explicit solution formulas for computing the intermediate filtering results we obtain an efficient algorithm. Inspired by [5] and [7], we propose the Multi-scale Graph-based Guided Filter (MG2F), which is to the best of our knowledge the first attempt of employing multi-scale graph representation as a guidance for an iterative graph spectral filtering in general and on CET data in particular. 1 Methodology We assume the noisy image Iη to be corrupted by white Gaussian noise, thus a suitable objective function would be Î f = arg min I f 1 2 ‖I f − Iη‖2, (1) ALBARQOUNI et al.: MULTI-SCALE GRAPH-BASED GUIDED FILTER (MG2F) 3 Figure 1: MG2F Framework: A noisy image slice from the 3D reconstructed tomogram is fed to the algorithm, where the graph is built on a selected scale space image (i.e. coarse grid) acting as a guidance for the regularized graph spectral filter. but we will augment this energy by a novel multi-scale graph regularisation as described in the following. 1.1 Graph Representation Given a noisy image Iη , we collect N overlapping patches covering the whole image P ∈ R √ n× √ n, which can be seen as data points ν = {ν1,ν2, ...,νN} ∈ Rn×N lying on a manifold M embedded in Rn space such that ν = EIη , where E is an operator collecting patches and vectorize it, cf. Figure 1. The relation between the data points can be represented by a k-NN connected, undirected, and weighted graph G = {ν ,ε,ω}, where ν is the data points (patches), ε is the set of edges, and ω is the set of edge weights. 1.1.1 Weight Assignment Assigning weights to the edges which exhibit a low SNR such in Cryo-ET data is challenging, therefore, we recall the scale-space theory [15] to build a Gaussian pyramid IGσs = Gσs ∗ Iη such that the noise manifests itself at certain structure scale σs and the semantical image appears clearly as shown in Figure 1, then the weights of the data points can be easily assigned using a heat kernel as follows: Wi j = exp− ‖νi−ν j‖22,σs σ2 h , εi j ∈ k-NN, 0, else. (2) Where ‖·‖2,σ is the Euclidean distance between two vectors at scale σ , however, σh is controlling the affinity of the neighbouring data points. Further, we denote the diagonal degree matrix by D, where Dii := ∑ j Wi j. 1.1.2 Graph Guidance Regularization In this paper, we are interested to preserve the intrinsic structure of the data points ν in the spectral filtering phase. It is worth mentioning that ν will be collected from the iterated 4 ALBARQOUNI et al.: MULTI-SCALE GRAPH-BASED GUIDED FILTER (MG2F) filtered image I f . We recall the definition of the Laplacian Quadratic Form [20], which can be represented as follows: Sσs(I f ) = ∑ (i, j)∈ε Wi j‖νi−ν j‖2 = 1 2 Tr ( νLσs ν T ) . (3) This expression can be interpreted as a regularization term that minimizes the distance between data points guided by, we denote it, the penalty Laplacian graph Lσs := D−W , which computed at different scales σs. The normalised Laplacian graph can be computed by ̃ Lσs := D − 1 2 Lσs D − 2 . The reader is referred to [8] for more details on graph operators. 1.2 Graph Spectral Filtering The spectrum of the graph σ(G) can be obtained from the eigenvalue decomposition of the normalised graph Laplacian ̃ Lσs :=UΛUT , where the eigenvalues Λ= diag{λ1,λ2, ...,λN}∈ [0,1] carry a notion of the graph frequencies, and the eigenvectors UT := {u1,u2, ...,uN} ∈ RN×n act as the orthogonal basis of the Graph Fourier Transform (GFT) [8], so we can write the transformed signal as follows ν̂ =UT ν . 1.2.1 Regularized Energy We define our objective function as follows: Î f = arg min I f { 1 2 ‖I f − Iη‖2 +αSσs(I f ) } , (4) where α > 0 is the regularization parameter. The solution can be written in a closed form: Î f = ET ( N ∑ i=1 1 (1+αλi) uiν̂i ) = ET ( 1 I +αL̃σs ) EIη , (5) where ET denotes the reshaping process of the previously vectorised patches. It becomes apparent from (5) that the signal is filtered on the spectral domain before doing the inverse GFT, where the spectral response of the filter h1(λi) = 1/(1+αλi) controls the frequency decay and thus the degree of smoothness. 1.3 Connection to Classical Filters Different classical filters can be expressed similar to (5) with different spectral filters, for instance, the Bilateral filter (BF) kernel can be written as νBF = D−1Wν , where its spectral response can be recast as a linear spectral filter h(λi) = (1− λi) [7], the same applies for non-local means filter (NLM), while the nonlinear anisotropic diffusion (NAD) has an exponential spectral filter h2(λi) = e−αλi . For the sake of having the power of diffusivity (fast decay) along with the regularized graph, we propose a new spectral filter h3(λi) = ei/(1+αλi), (6) where κ is a decaying factor. The proposed filter can be interpreted in the context of fractional derivative orders of Laplacian in Sobolev space [3], that shows a promising performance. A comparison of different filter responses is shown in Figure 2. ALBARQOUNI et al.: MULTI-SCALE GRAPH-BASED GUIDED FILTER (MG2F) 5 (a) h(λi) (b) h2(λi) (c) h3(λi) (d) Filter Response Figure 2: Spectral filters responses: (a) linear (i.e. Bilateral), (b) regularised graph, and (c) designed one, against the parameters λ (spectral frequency) and α (regularization parameter), (d) shows the line profile (α = 1) for different filters. 1.4 Stopping Criterion One can simply raise a question, why we need a stopping criterion where we have already a closed form solution for (4). Indeed, this optimal solution is designated for a specific scale, and since we are interested in having a multi-scale reconstruction, the resultant filtered image from the previous scale used as a guidance for the next scale, hence the need of a stopping criterion. Choosing it automatically is an important feature for variational approaches in general, [10] suggested one stopping criterion for Manifold de-noising, based on graph diffusion, therefore we employ the graph diffusion distance proposed by [9], ξ (Lk,Lk−1) := ‖e−L − e−Lk−1‖2F , which computes each iteration the distance between consequent graph Laplacians, which reflects the significant change in the filtering process. Then the optimization problem formulated as follows: Î f = arg min I f ,σs { 1 2 ‖I f − Iη‖2 +αSσs(I f ) } , s.t. ξ (Lk,Lk−1)≤ β , (7) where β is the desired distance and k is the iteration index, which can be minimized by Algorithm 1. Require: The noisy image Iη , patch size √ n, σs, k-NN, kmax, α and β . Ensure: The filtered image Î f in (7). 1: Initialize the Laplacian L0 = ones(N,N), If ← Iη 2: while ξ (Lk,Lk−1)> β or k < kmax do 3: Find the scale-space image IGσs = Gσs ∗ I k f . 4: σs← next finer level in the pyramid. 5: Collect patches (data points) ν = EIf . 6: Build the graph & assign weights for the k-NN patches νσs = EIGσs using (2). 7: Compute the normalised graph Laplacian ̃ Lσs and the graph spectrum σ(G). 8: Apply the spectral filtering Ik+1 f = E T (∑Ni=1 h(λi)uiν̂i). 9: Compute the graph diffusion ξ (Lk,Lk−1). 10: If ← I k+1 f 11: end while Algorithm 1: MG2F algorithm 6 ALBARQOUNI et al.: MULTI-SCALE GRAPH-BASED GUIDED FILTER (MG2F) (a) G.T. (b) Noisy (c) BF (d) BTR (e) EED (f) RGF (g) NAD (h) NLM (i) MG2F Algorithm PSNR MSE Parameters (dB) (10−3) Bilateral (BF) [19] (σi=0.5, σr=1.5, W=10) 17.49 174 Beltrami (BTR) [4] (δ=0.1, iter=10) 17.37 176 EED [18] (ρ=4, iter=30) 11.27 324 NAD [15] (iter=10, κ=0.3) 16.50 192 NLM [2] (P=7, W=21, σs=4σn ) 12.11 298 RGF [5] (σi=0.5, σr=1.5, iter=10) 17.49 174 MG2F (α=0.8, iter=4, σh=0.1) 17.78 169 Figure 3: Photographic Image: Results of different algorithms on Lena image(128X128, SNR=7) along with a tabulated comparison to the proposed MG2F filter. 2 Experiments and Results Our experiments are conducted on computer vision, simulated data, that rather mimics the complicated noise model in CET, as well as real CET tomographic data. We compare the results of our algorithm (MG2F) against common de-noising gold-standard filters in computer vision community, then we compare it with the successful filters in CET [12] such as Nonlinear Anisotropic Diffusion (NAD) [15] and Fast Non-Local Means (NLM) [2], and further with the recent scale-aware filter, the Rolling Guidance Filter (RGF) [5]. The filter parameters are tuned in an optimal way; either from the cited references or determined experimentally in order to visualise the feature of interest. Results are validated by different metrics; i.e. data with ground truth are validated using Peak Signal to Noise Ratio (PSNR) and Mean Square Error (MSE), however, we followed [6] and [12] for evaluating denoising methods on real data. Computer Vision: To give a good illustrative example, we run the algorithm on Lena image, which corrupted by an (i.i.d) Gaussian noise resulting in SNR of 7. Different algorithms are applied on this image, results are shown in Figure 3 for the cropped images. It is clear that our method gives an outperforming PSNR indicating for better contrast. Simulated Data: A GroEL tomogram, which is obtained from the Electron Microscopy Data Bank (EMDB-ID: 1AON), is generated using the TOM toolbox [13], where both (i.i.d) Gaussian and Poisson noise are added on the projection images resulting in a signal-to-noise ratio (SNR) of 0.1, then a slice of the reconstructed tomogram (the noisy image) is passed to different denoising algorithms. We run the algorithms on 150 random slices collected from 15 different tomograms. Results are validated by PSNR, our method shows significant performance (p < 0.01 by t-test, and p < 0.05 by Kolmogorov-Smirnov test), making the results consistently better compared to the other methods as shown in Figure 4(c). ALBARQOUNI et al.: MULTI-SCALE GRAPH-BASED GUIDED FILTER (MG2F) 7 (a) 64×64 (b) 128×128 (c) PSNR (d) FSC Figure 4: Sensitivity analysis: PSNR contour against k-NN and Patch size for different image sizes (a) 642 and (b) 1282. In (c) PSNR of different denoising algorithms (NLM, NAD, RGF and Ours respectively) for 150 slices (SNR=0.1) from 15 simulated tomograms. (d) FSC curve for different denoised 3D tomograms. (a) Noisy Image (b) NLM (c) NAD (d) RGF (e) MG2F Figure 5: 2D CET data: Filtering results on the tomogram along with the corresponding CNR of b) NLM (0.1979), c) NAD (0.2570), d) RGF (0.3146), e) Proposed MG2F (0.3150), where the arrows point to the fine structures on the membrane and the ellipse contains the inner core of HIV virus. Real CET Data: We also denoised an unstained CET HIV-1 virion (EMDB-ID: 1155), which can be considered as a benchmark data for de-noising in CET. Results are validated by the common validation measure in CET community for 2D images, namely, the Contrast to Noise ratio (CNR) [6]. Further, in the qualitative evaluation session with our clinical partners, they appreciated the enhanced contrast in our method, because the background in Figure 5(e) was carefully smoothed, while spiky signals such as membrane proteins (arrows) as well as the inner core of the HIV virus (ellipse) were well preserved. As the algorithm takes factors, such as patch size, redundancy, and scale-space, into account, we may conclude, based on their feedback and the resultant CNRs, that it is well suited for handling CET data. 3D Extension: Showing 3D data augments the computed statistics visually and gives a better interpretation for scientists, therefore, we extended the algorithm to handle 3D objects, which conceptually similar to our basic algorithm, however, blocks are collected from volumes rather than patches. An unstained HIV-1 tomogram is denoised using NAD (com8 ALBARQOUNI et al.: MULTI-SCALE GRAPH-BASED GUIDED FILTER (MG2F) (a) Noisy (b) BM3D (c) NAD (d) MG2F Figure 6: 3D CET Data: A comparison between different 3D filtering methods to our proposed MG2F method on real unstained HIV-1 data (EMDB-ID: 1155). monly used in CET), BM3D [1] (commonly used in Computer Vision) and our algorithm MG2F . The results are validated using both KL-divergence test (p < 0.1) and Fourier Shell Correlation (FSC) used in [12] as shown in Figure 6. The FSC curve shows the correlation of the corresponding frequency shells between the unprocessed/noisy and denoised tomograms, the blue line shows the auto-correlation of noisy data over the frequency, NAD has a smooth curve as expected due to the diffusion effect, BM3D is slightly better than MG2F in low frequencies, however, MG2F has a sharp decay in higher frequencies which is not the case for BM3D as shown in Figure 4(d). It is worth mentioning that the higher 0.5-cut-off frequency the higher resolution you get for these tomograms. Therefore, we can say MG2F performs better than NAD (has lower resolution), and BM3D (pass the higher frequencies). Sensitivity Analysis: Cross validation in general is a daunting task, however, it becomes more difficult when the feature assessment depends mainly on the experts opinion. Therefore, we performed a sensitivity analysis to investigate the effect of these hyperparamaters and suggest to update the scale-space parameter σs iteratively from the coarse level (given) to the fine level (until saturated). σh should be set based on the variance of data points, however, for the sake of simplicity, we normalize the data points before computing the weights. The impact of selecting different k-NN and Patch sizes (i.e. 3,5,7,9 and 11) on PSNR is shown in Figure 4. We observe that the algorithm converges at 3-7 iterations demonstrating the performance of the stopping criterion.